#### Start of replication of cross-CSO analyses ####


#*########################################*#
#### Data Preparation and descriptives #####
#*########################################*#

# Load packages
library(tidyverse)
library(lfe)
library(plm)
library(dplyr)
library(stargazer)
library(readxl)

# Clean environment
rm(list = ls())

# @USER: PLEASE CHANGE TO YOUR WORKING DIRECTORY
#setwd("C:/Users/HSmidt/Dropbox/Projekte/Paper_05 - CSO repression and Human Rights complaints/03_Code_and_Data/")
setwd("C:/Users/chrstein/Dropbox/CSO repression and Human Rights complaints/03_Code_and_Data/")

# Load data
load("10_replication_files/data_for_analysis_csos.RData")
ldata8d <-ldata8d %>% arrange(cowc, cso_name, year) %>% as.data.frame()
ldata <- ldata8d

# Check all CSOs that participated
unique(ldata$cso_name)

# Number survey responses by country
cso_dat <- ldata[!duplicated(ldata$cso_name), ]
tab <- table(cso_dat$cowc)
barplot(tab)

# Position of respondents
sort(table(cso_dat$position))

# Relation to government changed as consequence of complaint
table(cso_dat$relation_gov_changed)
categories <- c("Much worsened", "Worsened", "Not changed", "Improved",
                       "Much improved")
values <- c(5, 7, 27, 1, 0)
tab <- cbind(categories,values)
tab <- as.data.frame(tab)
tab$values <- as.numeric(tab$values)
category_order <- c("Much worsened", "Worsened", "Not changed", "Improved",
                    "Much improved") 
tab$categories <- factor(tab$categories, levels = category_order)


#### Figure 4 (main ms.) #####
ragg::agg_png("./10_replication_files/Outputs/Figures/Fig4_mainms_relationgov_changed.png"
               , width = 297, height = 210, units = "mm", res = 600)
ggplot(tab, aes(x = categories, y = values)) +
  theme_bw() + 
  theme(axis.text=element_text(size=16)
        ,axis.text.x=element_text(size=16, angle = 35, hjust = 1)
        ,panel.border = element_blank()) + ylim(c(0,30)) +
  geom_bar(stat = "identity", position="dodge", colour="black", alpha=0.2 ) +  xlab("") + ylab("") 
dev.off()


# Behavior of government changed as a consequence of complaint
table(cso_dat$behavior_gov_changed)
categories <- c("Much more repressive", "Somewhat more repressive", "Stayed the same", "Somewhat less repressive",
                "Much less repressive")
values <- c(3, 13, 24, 0, 0)
tab <- cbind(categories,values)
tab <- as.data.frame(tab)
tab$values <- as.numeric(tab$values)
category_order <-  c("Much more repressive", "Somewhat more repressive", "Stayed the same", "Somewhat less repressive",
                     "Much less repressive")
tab$categories <- factor(tab$categories, levels = category_order)

#### Figure 5 (main ms.) #####
ragg::agg_png("./10_replication_files/Outputs/Figures/Fig5_mainms_behaviorgov_changed.png"
               , width = 297, height = 210, units = "mm", res = 600)
ggplot(tab, aes(x = categories, y = values)) +
  theme_bw() + 
  theme(plot.margin = margin(1,1,0.2,1, "cm")
        ,axis.text.x=element_text(size=16, angle = 35, hjust = 1)
        ,axis.title=element_text(size=16)
        ,axis.text.y = element_text(size=16)
        ,panel.border = element_blank()) + 
  geom_bar(stat = "identity", position="dodge", colour="black", alpha=0.2) +  xlab("") + ylab("") 
dev.off()

# Distribution of main dependent variable
table(ldata$relation_gov_yearly)
table(ldata$relation_gov_yearly_num)
categories <- c("Very negative", "Negative", "Neutral", "Positive", "Very positive")
values <- as.data.frame(table(ldata$relation_gov_yearly_num))[,2]
tab <- cbind(categories,values)
tab <- as.data.frame(tab)
tab$values <- as.numeric(tab$values)
category_order <-  c("Very negative", "Negative", "Neutral", "Positive", "Very positive")
tab$categories <- factor(tab$categories, levels = category_order)

#### Figure A.3 (appendix) #####
ragg::agg_png("./10_replication_files/Outputs/Figures/FigA3_AppD_behaviorgov_changed.png"
               , width = 297, height = 210, units = "mm", res = 600)
#pdf(file = "./04_analysis_cross_csos/outputs/relation_to_gov.pdf")
ggplot(tab, aes(x = categories, y = values)) +
  theme_bw() + 
  theme(plot.margin = margin(1,1,1,1, "cm")
        ,axis.text.y = element_text(size=16)
        ,axis.text.x=element_text(size=16, angle = 0, hjust = 0.5)
        ,axis.title=element_text(size=16)
        , axis.title.y = element_text(size=16)
        ,panel.border = element_blank()) + 
  geom_bar(stat = "identity", position="dodge", colour="black", alpha=0.2) + xlab("") + ylab("Number of CSO-years") 
dev.off()


# Known that it has been mentioned in Communication
table(cso_dat$mentioned_known)

# Known that it has filed complaint
table(cso_dat$filed_known)

 

#*###############################################*#
#### Create relevant variables for the analyses ####
#*###############################################*#
 
# Years_filed_long = 1 in years where CSO i filed complaint
ldata$years_filed_long <- ifelse(ldata$years_filed_long>0,1,0)

# Create binary variable "treatment" that is 1 in first year of communication and all subsequent years
# Create binary variable "complaint" that is 1 in first year of complaint and all subsequent years
# Create count variable "complaint_years_since" that counts years since first complaint year + 1
# Create count variable "treatment_years_since" that counts years since first communication year + 1
ldata <- ldata %>% arrange( cso_name, year ) %>% group_by(cso_name) %>% 
                                mutate(treatment = cumsum(year_firstcomm),
                                       complaint = cumsum(years_filed_long),
                                       complaint = ifelse(complaint>0,1,0),
                                       complaint_years_since = cumsum(complaint),
                                       treatment_years_since = cumsum(treatment)) 


## Make lag
ldata <- ldata %>% group_by(cowc, cso_name) %>%
              arrange(cowc, cso_name, year) %>%
              mutate( relation_gov_yearly_num_lag1 = dplyr::lag(relation_gov_yearly_num,1),
                      year_firstcomm_lag1 = dplyr::lag(year_firstcomm,1),
                      treatment_lag1 = dplyr::lag(treatment, 1),
                      complaint_lag1 = dplyr::lag(complaint, 1),
                      treatment_years_since_lag1 = dplyr::lag(treatment_years_since, 1),
                      complaint_years_since_lag1 = dplyr::lag(complaint_years_since, 1),
                      aid_odapercap_lag = dplyr::lag(aid_odapercap,1),
                      new_entries_sum_lag1 = dplyr::lag(new_entries_sum, 1),
                      new_entries_sum_lag2 = dplyr::lag(new_entries_sum, 2),
                      Rapporteurs_active_number_lag1 = dplyr::lag(Rapporteurs_active_number, 1),
                      Rapporteurs_active_number_lag2 = dplyr::lag(Rapporteurs_active_number, 2),
                      UNSP_language_n_lag1 = dplyr::lag(UNSP_language_n,1),
                      UNSP_language_n_lag2 = dplyr::lag(UNSP_language_n,2),
                      spvisits_lag1 = dplyr::lag(spvisits,1) )


# Set period of analyses to 2010-2022
ldata <- subset(ldata, year>=2010 & year<=2022)

# Complete Cases?
complete_ldata <- ldata[complete.cases(ldata[c("treatment","spvisits_lag1","new_entries_sum","Rapporteurs_active_number","v2peedueq_lag1","conflict_lag","v2mecenefm_lag1","v2juhcind_lag1","log_gdpcap_lag1")]),]

# How many years, countries, and CSOs in sample
nrow(complete_ldata)
table(ldata$year)
table(ldata$cowc)
unique(ldata$cso_name)

###################################
#### Models  in Main body of Ms. ##
###################################

controls <- paste0("+ v2peedueq_lag1 + conflict_lag + v2mecenefm_lag1 + v2juhcind_lag1 + log_gdpcap_lag1")

# FE models

ldata$years_since <- ldata$treatment_years_since

# Binary variable for year of communication and all subsequent years
f0 <- as.formula( paste0("relation_gov_yearly_num  ~ treatment | 0 | 0 | cso_name + year"))
f1 <- as.formula( paste0("relation_gov_yearly_num ~  treatment | cso_name | 0 | cso_name + year"))
f2 <- as.formula( paste0("relation_gov_yearly_num ~  treatment ", controls,  " + years_since | cso_name | 0 | cso_name + year"))
f3 <- as.formula( paste0("relation_gov_yearly_num ~  treatment ", controls, " + years_since + lag(relation_gov_yearly_num,1) | lag(relation_gov_yearly_num,2:99)") )
f4 <- as.formula( paste0("relation_gov_yearly_num ~  " , controls, " + years_since + spvisits_lag1 | cso_name | (treatment ~  new_entries_sum + Rapporteurs_active_number) | cso_name + year"))

m0 <- felm(formula = f0, data = ldata )
m1 <- felm(formula = f1, data = ldata )
m2 <- felm(formula = f2, data = ldata )
m3 <- pgmm(formula = f3, data = ldata, index=c("cso_name", "year"), transformation = 'ld',  model = 'twostep', effect = 'individual')
m4 <- felm(formula = f4, data = ldata )

#*###########################################*#
#### First Stage Model in Main Body of Ms. ####
#*###########################################*#
f4_fs <- as.formula(paste0("treatment ~ new_entries_sum + Rapporteurs_active_number", controls,  "+ spvisits_lag1 + years_since | cso_name | 0 | cso_name"))
m4_fs <- felm(formula = f4_fs, data = model.frame(m4))
summary(m4_fs)
#save(m4_fs, file="04_analysis_cross_csos/first_stage_models_csos.R")

## Get fstat
fstat2 <- round(summary(m4_fs)$fstat, 3)
fstat2

# Only first stage on CSO-level
custom_significance <- c("\\dag", "*", "**", "***")
star_first_full <- stargazer(m4_fs
                             , title = "First stage model"
                             , digits=3
                             , font.size="scriptsize"
                             , omit.stat = c("n", "ser")
                             , add.lines = list(c("Number of observations"
                                                  , length(m4_fs$fitted.values)) ,
                                                c("Country fixed effects", "Yes"),
                                                c("Time-varying controls", "Yes"),
                                                c("Lagged depend. var.", "No"),
                                                c("First stage F-statistic", fstat2))
                             , omit = c("v2peedueq_lag1",
                                        "conflict_lag",
                                        "v2mecenefm_lag1",
                                        "v2juhcind_lag1",
                                        "log_gdpcap_lag1",
                                        "spvisits_lag1",
                                        "years_since")
                             , covariate.labels = c( "New UNSP experts"
                                                     , "Number of UNSP experts")
                             , column.labels = c("(1)")
                             , dep.var.labels.include = F
                             , dep.var.caption = "Dependent variable: UNSP communications"
                             , column.sep.width = "1pt"
                             , align=F
                             , label = "tab:first_stage_CSO"
                             , model.numbers=F, model.names=F
                             , table.placement="H"
                             , star.char = custom_significance
                             , star.cutoffs = c(0.1, 0.05, 0.01, 0.001))

note.latex <- "\\multicolumn{2}{l} {\\parbox[t]{9.5cm}{ \\textit{Notes:} The models are estimated with ordinary least squares regression. Cluster-robust standard errors are in parentheses. $^{\\dag}$p$<$0.1; $^{*}$p$<$0.05; $^{**}$p$<$0.01; $^{***}$p$<$0.001}.} \\\\"
star_first_full[grepl("Note",star_first_full)] <- note.latex
cat(star_first_full, sep = "\n")
#write(star_first_full, "./02_analysis_cross_national/outputs/first_stage_CSO.tex")
write(star_first_full, "./10_replication_files/Outputs/Tables/Table2_mainms_first_stage_CSO.tex")

# Summary stats
fstat <- round(summary(m4_fs)$fstat, 3)
fstat
sargan <- plm::sargan(m3)
sargan
autocorr1 <- mtest(m3,1)
autocorr1
autocorr2 <- mtest(m3,2)
autocorr2  
autocorr3 <- mtest(m3,3)
autocorr3  

# Main body of paper
custom_significance <- c("\\dag", "*", "**", "***")
star_reduced <- stargazer(m0, m1, m2, m3, m4
                          , type="latex"
                          , title = "Reported CSO-government relations and UNSP communications"
                          , digits=3
                          , font.size="scriptsize"
                          , notes.align ="l"
                          , omit.stat = c("n", "ser", "rsq", "adj.rsq")
                          , add.lines = list(c("Number of observations"
                                               , length(m0$fitted.values)
                                               , length(m1$fitted.values)
                                               , length(m2$fitted.values)
                                               , length(m2$fitted.values)
                                               , length(m4$fitted.values)) ,
                                             c("CSO fixed effects", "No", "Yes", "Yes", "Yes", "Yes"),
                                             c("Time-varying controls", "No", "No", "Yes", "Yes", "Yes"),
                                             c("Lagged depend. var.", "No", "No", "No", "Yes", "No"),
                                             c("First stage F-statistic", "NA", "NA", "NA", "NA",  fstat))
                          , keep =  c("treatment")
                          , covariate.labels = c("Post-communication"
                                                 ,"Post-communication (instr.)")
                          , column.labels = c("(1)","(2)","(3)","(4)","(5)")
                          , dep.var.labels.include = F
                          , dep.var.caption = "DV: Self-reported relation with government (original survey)"
                          , column.sep.width = "1pt"
                          , align=F
                          , model.numbers=F, model.names=F
                          , label="tab:relation_gov_yearly_reduced"
                          , table.placement="H"
                          , star.char = custom_significance
                          , star.cutoffs = c(0.1, 0.05, 0.01, 0.001))


note.latex <- "\\multicolumn{6}{l} {\\parbox[t]{12cm}{ \\textit{Notes:} The models in columns 1, 2, 3, and 5 are estimated with ordinary least squares regression, while the model in column 4 is estimated using generalized methods of moments. The lagged dependent variable in the model in column 4 is instrumented with deeper lags. In the models in columns 3, 4 and 5, we control for education, conflict, media freedom, judicial independence, and logged GDP per capita, all lagged one year, as well as number of years since the first communication targeting the CSO. Cluster-robust standard errors are in parentheses. $^{\\dag}$p$<$0.1; $^{*}$p$<$0.05; $^{**}$p$<$0.01; $^{***}$p$<$0.001}.} \\\\"
star_reduced[grepl("Note",star_reduced)] <- note.latex
cat(star_reduced, sep = "\n")
#write(star_reduced, "./04_analysis_cross_csos/outputs/relation_gov_yearly_reduced.tex")

#### Table A.18 for Appendix D ####
star_full <- stargazer(m0, m1, m2, m3, m4
                          , type="latex"
                          , title = "Reported CSO-government relations and UNSP communications"
                          , digits=3
                          , font.size="scriptsize"
                          , notes.align ="l"
                          , omit.stat = c("n", "ser", "rsq", "adj.rsq")
                          , add.lines = list(c("Number of observations"
                                               , length(m0$fitted.values)
                                               , length(m1$fitted.values)
                                               , length(m2$fitted.values)
                                               , length(m2$fitted.values)
                                               , length(m4$fitted.values)) ,
                                             c("CSO fixed effects", "No", "Yes", "Yes", "Yes", "Yes"),
                                             c("Time-varying controls", "No", "No", "Yes", "Yes", "Yes"),
                                             c("Lagged depend. var.", "No", "No", "No", "Yes", "No"),
                                             c("First stage F-statistic", "NA", "NA", "NA", "NA",  fstat))
                           , covariate.labels = c("Post-communication years"
                                                 , "Education (lag 1 yr)"
                                                 , "Conflict (lag 1 yr)"
                                                 , "Freedom of media (lag 1 yr)"
                                                 , "Judicial independence (lag 1 yr)"
                                                 , "Logged GDP per capita (lag 1 yr)"
                                                 , "Years since first communication"
                                                 , "Relation with government (lag 1 yr, instr.)"
                                                 , "UNSP country visits"
                                                 , "Post-communication (instr.)"
                                                 , "Constant")
                          , column.labels = c("(1)","(2)","(3)","(4)","(5)")
                          , dep.var.labels.include = F
                          , dep.var.caption = "DV: Self-reported relation with government (original survey)"
                          , column.sep.width = "1pt"
                          , align=F
                          , model.numbers=F, model.names=F
                          , label="tab:relation_gov_yearly_full"
                          , table.placement="H"
                          , star.char = custom_significance
                          , star.cutoffs = c(0.1, 0.05, 0.01, 0.001))


note.latex <- "\\multicolumn{6}{l} {\\parbox[t]{14cm}{ \\textit{Notes:} The models in columns 1, 2, 3, and 5 are estimated with ordinary least squares regression, while the model in column 4 is estimated using generalized methods of moments. The lagged dependent variable in the model in column 4 is instrumented with deeper lags. Cluster-robust standard errors are in parentheses. $^{\\dag}$p$<$0.1; $^{*}$p$<$0.05; $^{**}$p$<$0.01; $^{***}$p$<$0.001}.} \\\\"
star_full[grepl("Note",star_full)] <- note.latex
cat(star_full, sep = "\n")
#write(star_full, "./04_analysis_cross_csos/outputs/relation_gov_yearly_full.tex")
write(star_full, "./10_replication_files/Outputs/Tables/TableA18_AppD_relation_gov_yearly_full.tex")



## Substantive interpretation

# Effects
effect_sizes_csos <- matrix(NA, nrow=5, ncol=5)
coef_m0 <- MASS::mvrnorm(1000, coefficients(m0), vcov(m0))[,2] / sd( model.frame(m0)[,c("relation_gov_yearly_num")] )
effect_sizes_csos[1,] <- quantile(coef_m0, c(0.025, 0.05, 0.5, 0.95, 0.975))
coef_m1 <- MASS::mvrnorm(1000, coefficients(m1), vcov(m1))[,1] / sd( model.frame(m1)[,c("relation_gov_yearly_num")] )
effect_sizes_csos[2,] <- quantile(coef_m1, c(0.025, 0.05, 0.5, 0.95, 0.975))
coef_m2 <- MASS::mvrnorm(1000, coefficients(m2), vcov(m2))[,1] / sd( model.frame(m2)[,c("relation_gov_yearly_num")] )
effect_sizes_csos[3,] <- quantile(coef_m2, c(0.025, 0.05, 0.5, 0.95, 0.975))
coef_m3 <- MASS::mvrnorm(1000, coefficients(m3), vcov(m3))[,1] / sd( model.frame(m2)[,c("relation_gov_yearly_num")] )
effect_sizes_csos[4,] <- quantile(coef_m3, c(0.025, 0.05, 0.5, 0.95, 0.975))
coef_m4 <- MASS::mvrnorm(1000, coefficients(m4), vcov(m4))[,8] / sd( model.frame(m4)[,c("relation_gov_yearly_num")] )
effect_sizes_csos[5,] <- quantile(coef_m4, c(0.025, 0.05, 0.5, 0.95, 0.975))

# Compared to effect of armed conflict
summary(m4)
coefficients(m4)[2] / sd( model.frame(m4)[,c("relation_gov_yearly_num")] ) 



##### Figure 6 (main ms) #####
ragg::agg_png("./10_replication_files/Outputs/Figures/Fig6_mainms_marg_effect_plot_crosscso.png"
               , width = 297, height = 210, units = "mm", res = 600)
par(mar = c(5.5, 16, 2.1, 2.8), mgp = c(3.5, 1, 0), xpd = TRUE)
plot( y = c(1,2,3,4,5)
      , x = rev(effect_sizes_csos[,3])*100
      , pch= rev(c(16,15,17,18,20))
      , cex=1.5
      , type="p"
      , cex.lab = 1.5
      , axes= F, xlab = "Effect of post-UNSP communication period (as % st. dev. of CSO-gov. relations)"
      , ylab = ""
      , ylim = c(0.5, 5.5), xlim = c(-150,50) )
arrows( x0 = rev(effect_sizes_csos[,1]*100), x1 = rev(effect_sizes_csos[,5]*100)
        , y0 = c(1,2,3,4,5),      y1 = c(1,2,3,4,5)
        , angle=90, length=0, lwd=1)
arrows( x0 = rev(effect_sizes_csos[,2]*100), x1 = rev(effect_sizes_csos[,4]*100)
        , y0 = c(1,2,3,4,5),      y1 = c(1,2,3,4,5)
        , angle=90, length=0, lwd=2)
lines(y=c(0.5,5.5), x = c(0,0), lty="dashed")
axis(1, at = seq(-150, 50, by=50), labels = seq(-150, 50, by=50),  cex.axis=1.5)
text(x = rep(-240,5), y=c(1,2,3,4,5), labels=c("With CSO FEs, \ncontrols, and instrument"
                                               , "With CSO FEs, \ncontrols, and instrumented \nlagged DV"
                                               , "With CSO FEs and \ncontrols"
                                               , "With CSO FEs"
                                               , "Bivariate regression"), cex=1.5, pos = 4)
text(x= -240, y=5.5, labels = "Models:", font=2, cex=1.5, pos=4)
dev.off()


## Test sensitivity of IV
Y <- "relation_gov_yearly_num" # Y: outcome of interest
D <-"treatment" # D: endogenous treatment
Z <- c("new_entries_sum", "Rapporteurs_active_number" ) # Z: instrumental variable; add: "Rapporteurs_active_number"
controls <- c("v2peedueq_lag1", "conflict_lag", "v2mecenefm_lag1", "v2juhcind_lag1", "log_gdpcap_lag1", "spvisits_lag1", "treatment_years_since") # covariates of control variables
FE <- c("cso_name") # FE
cl <- c("cso_name") # clusters
 
library(ivDiag)
g <- ivDiag(data = ldata, Y=Y, D = D, Z = Z, controls = controls, FE= FE, cl = cl)
g$est_ols
g$est_2sls
g$AR
g$F_stat
g$rho
g$est_fs
g$est_rf


#### Re-run the analysis for subset of domestic NGOs ####

# Create dataset without the international NGO in the sample
ldata_domestic <- ldata %>% filter(cso_name != "Lawyers for Human Rights")

# Binary variable for year of communication and all subsequent years
m0 <- felm(formula = f0, data = ldata_domestic)
m1 <- felm(formula = f1, data = ldata_domestic)
m2 <- felm(formula = f2, data = ldata_domestic)
m3 <- pgmm(formula = f3, data = ldata_domestic, index=c("cso_name", "year"), transformation = 'ld',  model = 'twostep', effect = 'individual')
m4 <- felm(formula = f4, data = ldata_domestic)

# Summary stats
fstat <- round(summary(m4_fs)$fstat, 3)
fstat
sargan <- plm::sargan(m3)
sargan
autocorr1 <- mtest(m3,1)
autocorr1
autocorr2 <- mtest(m3,2)
autocorr2  
autocorr3 <- mtest(m3,3)
autocorr3  

# Main body of paper
custom_significance <- c("\\dag", "*", "**", "***")
star_reduced <- stargazer(m0, m1, m2, m3, m4
                          , type="latex"
                          , title = "Reported CSO-government relations and UNSP communications"
                          , digits=3
                          , font.size="scriptsize"
                          , notes.align ="l"
                          , omit.stat = c("n", "ser", "rsq", "adj.rsq")
                          , add.lines = list(c("Number of observations"
                                               , length(m0$fitted.values)
                                               , length(m1$fitted.values)
                                               , length(m2$fitted.values)
                                               , length(m2$fitted.values)
                                               , length(m4$fitted.values)) ,
                                             c("CSO fixed effects", "No", "Yes", "Yes", "Yes", "Yes"),
                                             c("Time-varying controls", "No", "No", "Yes", "Yes", "Yes"),
                                             c("Lagged depend. var.", "No", "No", "No", "Yes", "No"),
                                             c("First stage F-statistic", "NA", "NA", "NA", "NA",  fstat))
                          , keep =  c("treatment")
                          , covariate.labels = c("Post-communication"
                                                 ,"Post-communication (instr.)")
                          , column.labels = c("(1)","(2)","(3)","(4)","(5)")
                          , dep.var.labels.include = F
                          , dep.var.caption = "DV: Self-reported relation with government (original survey)"
                          , column.sep.width = "1pt"
                          , align=F
                          , model.numbers=F, model.names=F
                          , label="tab:relation_gov_yearly_reduced_onlydomestic"
                          , table.placement="H"
                          , star.char = custom_significance
                          , star.cutoffs = c(0.1, 0.05, 0.01, 0.001))


note.latex <- "\\multicolumn{6}{l} {\\parbox[t]{14cm}{ \\textit{Notes:} The models in columns 1, 2, 3, and 5 are estimated with ordinary least squares regression, while the model in column 4 is estimated using generalized methods of moments. The lagged dependent variable in the model in column 4 is instrumented with deeper lags. In the models in columns 3, 4 and 5, we control for education, conflict, media freedom, judicial independence, and logged GDP per capita, all lagged one year, as well as number of years since the first communication targeting the CSO. Cluster-robust standard errors are in parentheses. $^{\\dag}$p$<$0.1; $^{*}$p$<$0.05; $^{**}$p$<$0.01; $^{***}$p$<$0.001}.} \\\\"
star_reduced[grepl("Note",star_reduced)] <- note.latex
cat(star_reduced, sep = "\n")
#write(star_reduced, "./04_analysis_cross_csos/outputs/relation_gov_yearly_reduced_onlydomestic.tex")


# Appendix
star_full <- stargazer(m0, m1, m2, m3, m4
                       , type="latex"
                       , title = "Reported CSO-government relations and UNSP communications"
                       , digits=3
                       , font.size="scriptsize"
                       , notes.align ="l"
                       , omit.stat = c("n", "ser", "rsq", "adj.rsq")
                       , add.lines = list(c("Number of observations"
                                            , length(m0$fitted.values)
                                            , length(m1$fitted.values)
                                            , length(m2$fitted.values)
                                            , length(m2$fitted.values)
                                            , length(m4$fitted.values)) ,
                                          c("Country fixed effects", "No", "Yes", "Yes", "Yes", "Yes"),
                                          c("Time-varying controls", "No", "No", "Yes", "Yes", "Yes"),
                                          c("Lagged depend. var.", "No", "No", "No", "Yes", "No"),
                                          c("First stage F-statistic", "NA", "NA", "NA", "NA",  fstat))
                       , covariate.labels = c("Post-communication years"
                                              , "Education (lag 1 yr)"
                                              , "Conflict (lag 1 yr)"
                                              , "Freedom of media (lag 1 yr)"
                                              , "Judicial independence (lag 1 yr)"
                                              , "GDP per capita (lag 1 yr)"
                                              , "Years since first communication"
                                              , "CSO-reported repression (lag 1 yr, instr.)"
                                              , "UNSP country visits"
                                              , "Post-communication (instr.)")
                       , column.labels = c("(1)","(2)","(3)","(4)","(5)")
                       , dep.var.labels.include = F
                       , dep.var.caption = "DV: Self-reported relation with government (original survey)"
                       , column.sep.width = "1pt"
                       , align=F
                       , model.numbers=F, model.names=F
                       , label="tab:relation_gov_yearly_full_onlydomestic"
                       , table.placement="H"
                       , star.char = custom_significance
                       , star.cutoffs = c(0.1, 0.05, 0.01, 0.001))


note.latex <- "\\multicolumn{6}{l} {\\parbox[t]{14cm}{ \\textit{Notes:} The models in columns 1, 2, 3, and 5 are estimated with ordinary least squares regression, while the model in column 4 is estimated using generalized methods of moments. The lagged dependent variable in the model in column 4 is instrumented with deeper lags. Cluster-robust standard errors are in parentheses. $^{\\dag}$p$<$0.1; $^{*}$p$<$0.05; $^{**}$p$<$0.01; $^{***}$p$<$0.001}.} \\\\"
star_full[grepl("Note",star_full)] <- note.latex
cat(star_full, sep = "\n")
#write(star_full, "./04_analysis_cross_csos/outputs/relation_gov_yearly_full_onlydomestic.tex")


#### Between CSO analysis, for the time since Communications ####

# Load dataset that contains information on all mentioned years (not only first mentioned)
data_2 <- read_excel("10_replication_files/data_for_analysis_csos_mentioned_years.xlsx")

# Filter relevant variables and merge them to analysis dataset
data_2 <- select(data_2, cso_name, year, year_mentioned)
data_2$cso_year <- paste(data_2$cso_name, "_", data_2$year)
ldata$cso_year <- paste(ldata$cso_name, "_", ldata$year)
data_2 <- select(data_2, cso_year, year_mentioned)
ldata <- left_join(ldata, data_2, by = "cso_year")

# Add a counter that begins at 1 for the most recent communication for a given CSO
ldata <- ldata %>% 
  group_by(cso_name) %>% 
  mutate(lastyear = max(year[year_mentioned == 1]),
         since_last_mentioned = ifelse(year < lastyear,0, year - lastyear + 1 )) %>% 
  ungroup()

# Filter only the current ratings
onlymostrecent <- ldata %>% filter(year == 2022)

# Code CSO level covariates
onlymostrecent$human_rights_cso <- 0
onlymostrecent$human_rights_cso <- ifelse(onlymostrecent$cso_name == "Gulf Centre for Human Rights", 1, onlymostrecent$human_rights_cso)
onlymostrecent$human_rights_cso <- ifelse(onlymostrecent$cso_name == "Belarusian Helsinki Committee", 1, onlymostrecent$human_rights_cso)
onlymostrecent$human_rights_cso <- ifelse(onlymostrecent$cso_name == "Ukrainian Helsinki Human Rights Union", 1, onlymostrecent$human_rights_cso)
onlymostrecent$human_rights_cso <- ifelse(onlymostrecent$cso_name == "MADPET-Malaysians Against Death Penalty and Torture", 1, onlymostrecent$human_rights_cso)
onlymostrecent$human_rights_cso <- ifelse(onlymostrecent$cso_name == "Myanmar Ethnic Rohingya Human Rights Organization in Malaysia (MERHROM)", 1, onlymostrecent$human_rights_cso)
onlymostrecent$human_rights_cso <- ifelse(onlymostrecent$cso_name == "Pusat KOMAS", 1, onlymostrecent$human_rights_cso)
onlymostrecent$human_rights_cso <- ifelse(onlymostrecent$cso_name == "A tu encuentro guanajuato", 1, onlymostrecent$human_rights_cso)
onlymostrecent$human_rights_cso <- ifelse(onlymostrecent$cso_name == "Colectivo Buscadoras Guanajuato", 1, onlymostrecent$human_rights_cso)
onlymostrecent$human_rights_cso <- ifelse(onlymostrecent$cso_name == "CODESA", 1, onlymostrecent$human_rights_cso)
onlymostrecent$human_rights_cso <- ifelse(onlymostrecent$cso_name == "Legal Defence and Assistance Project (LEDAP)", 1, onlymostrecent$human_rights_cso)
onlymostrecent$human_rights_cso <- ifelse(onlymostrecent$cso_name == "Center for International Law (CenterLaw)", 1, onlymostrecent$human_rights_cso)
onlymostrecent$human_rights_cso <- ifelse(onlymostrecent$cso_name == "Karapatan Alliance Philippines", 1, onlymostrecent$human_rights_cso)
onlymostrecent$human_rights_cso <- ifelse(onlymostrecent$cso_name == "National Union of Peoples' Lawyers", 1, onlymostrecent$human_rights_cso)
onlymostrecent$human_rights_cso <- ifelse(onlymostrecent$cso_name == "Task Force Detainees of the Philippines (TFDP)", 1, onlymostrecent$human_rights_cso)
onlymostrecent$human_rights_cso <- ifelse(onlymostrecent$cso_name == "Tongtongan ti Umili", 1, onlymostrecent$human_rights_cso)
onlymostrecent$human_rights_cso <- ifelse(onlymostrecent$cso_name == "SOS Racismo", 1, onlymostrecent$human_rights_cso)
onlymostrecent$human_rights_cso <- ifelse(onlymostrecent$cso_name == "Òmnium Cultural", 1, onlymostrecent$human_rights_cso)
onlymostrecent$human_rights_cso <- ifelse(onlymostrecent$cso_name == "INFORM Human Rights Documentation centre", 1, onlymostrecent$human_rights_cso)
onlymostrecent$human_rights_cso <- ifelse(onlymostrecent$cso_name == "Lawyers for Human Rights", 1, onlymostrecent$human_rights_cso)
onlymostrecent$human_rights_cso <- ifelse(onlymostrecent$cso_name == "Right to Life Human Rights Centre", 1, onlymostrecent$human_rights_cso)
onlymostrecent$human_rights_cso <- ifelse(onlymostrecent$cso_name == "Association for Monitoring Equal Rights", 1, onlymostrecent$human_rights_cso)
onlymostrecent$human_rights_cso <- ifelse(onlymostrecent$cso_name == "DefendDefenders (East and Horn of Africa Human Rights Defenders Project)", 1, onlymostrecent$human_rights_cso)
onlymostrecent$human_rights_cso <- ifelse(onlymostrecent$cso_name == "Nakivale Victims Association", 1, onlymostrecent$human_rights_cso)
onlymostrecent$english_website <- 0
onlymostrecent$english_website <- ifelse(onlymostrecent$cso_name == "Gulf Centre for Human Rights", 1, onlymostrecent$english_website)
onlymostrecent$english_website <- ifelse(onlymostrecent$cso_name == "Belarusian Helsinki Committee", 1, onlymostrecent$english_website)
onlymostrecent$english_website <- ifelse(onlymostrecent$cso_name == "NGO Ecohome", 1, onlymostrecent$english_website)
onlymostrecent$english_website <- ifelse(onlymostrecent$cso_name == "Ukrainian Helsinki Human Rights Union", 1, onlymostrecent$english_website)
onlymostrecent$english_website <- ifelse(onlymostrecent$cso_name == "ZVYANO", 1, onlymostrecent$english_website)
onlymostrecent$english_website <- ifelse(onlymostrecent$cso_name == "All IndiaUnion of Forest Working People", 1, onlymostrecent$english_website)
onlymostrecent$english_website <- ifelse(onlymostrecent$cso_name == "Banglar Manabadhikar Suraksha Mancha (MASUM)", 1, onlymostrecent$english_website)
onlymostrecent$english_website <- ifelse(onlymostrecent$cso_name == "InformAction Kenya", 1, onlymostrecent$english_website)
onlymostrecent$english_website <- ifelse(onlymostrecent$cso_name == "Strategies for Northern Development", 1, onlymostrecent$english_website)
onlymostrecent$english_website <- ifelse(onlymostrecent$cso_name == "Malaysian Trade Union Congress", 1, onlymostrecent$english_website)
onlymostrecent$english_website <- ifelse(onlymostrecent$cso_name == "Myanmar Ethnic Rohingya Human Rights Organization in Malaysia (MERHROM)", 1, onlymostrecent$english_website)
onlymostrecent$english_website <- ifelse(onlymostrecent$cso_name == "Persatuan Sahabat Wanita Selangor", 1, onlymostrecent$english_website)
onlymostrecent$english_website <- ifelse(onlymostrecent$cso_name == "Pertubuhan Pelindung Khazanah Alam", 1, onlymostrecent$english_website)
onlymostrecent$english_website <- ifelse(onlymostrecent$cso_name == "Pusat KOMAS", 1, onlymostrecent$english_website)
onlymostrecent$english_website <- ifelse(onlymostrecent$cso_name == "SIS Forum Malaysia", 1, onlymostrecent$english_website)
onlymostrecent$english_website <- ifelse(onlymostrecent$cso_name == "Sabah Women's Action-Resource Group or SAWO", 1, onlymostrecent$english_website)
onlymostrecent$english_website <- ifelse(onlymostrecent$cso_name == "The Christian Federation of Malaysia", 1, onlymostrecent$english_website)
onlymostrecent$english_website <- ifelse(onlymostrecent$cso_name == "Amazigh Network for Citizenship - Azta Amazigh", 1, onlymostrecent$english_website)
onlymostrecent$english_website <- ifelse(onlymostrecent$cso_name == "Legal Defence and Assistance Project (LEDAP)", 1, onlymostrecent$english_website)
onlymostrecent$english_website <- ifelse(onlymostrecent$cso_name == "AGHAM Advocates of Science and Technology for the People", 1, onlymostrecent$english_website)
onlymostrecent$english_website <- ifelse(onlymostrecent$cso_name == "Center for International Law (CenterLaw)", 1, onlymostrecent$english_website)
onlymostrecent$english_website <- ifelse(onlymostrecent$cso_name == "GABRIELA", 1, onlymostrecent$english_website)
onlymostrecent$english_website <- ifelse(onlymostrecent$cso_name == "Karapatan Alliance Philippines", 1, onlymostrecent$english_website)
onlymostrecent$english_website <- ifelse(onlymostrecent$cso_name == "National Union of Peoples' Lawyers", 1, onlymostrecent$english_website)
onlymostrecent$english_website <- ifelse(onlymostrecent$cso_name == "Task Force Detainees of the Philippines (TFDP)", 1, onlymostrecent$english_website)
onlymostrecent$english_website <- ifelse(onlymostrecent$cso_name == "Tongtongan ti Umili", 1, onlymostrecent$english_website)
onlymostrecent$english_website <- ifelse(onlymostrecent$cso_name == "Fundación Secretariado Gitano", 1, onlymostrecent$english_website)
onlymostrecent$english_website <- ifelse(onlymostrecent$cso_name == "The Catalan National Assembly", 1, onlymostrecent$english_website)
onlymostrecent$english_website <- ifelse(onlymostrecent$cso_name == "Òmnium Cultural", 1, onlymostrecent$english_website)
onlymostrecent$english_website <- ifelse(onlymostrecent$cso_name == "Center for Women & Development", 1, onlymostrecent$english_website)
onlymostrecent$english_website <- ifelse(onlymostrecent$cso_name == "INFORM Human Rights Documentation centre", 1, onlymostrecent$english_website)
onlymostrecent$english_website <- ifelse(onlymostrecent$cso_name == "Lawyers for Human Rights", 1, onlymostrecent$english_website)
onlymostrecent$english_website <- ifelse(onlymostrecent$cso_name == "Pupil Salvation Forum", 1, onlymostrecent$english_website)
onlymostrecent$english_website <- ifelse(onlymostrecent$cso_name == "Right to Life Human Rights Centre", 1, onlymostrecent$english_website)
onlymostrecent$english_website <- ifelse(onlymostrecent$cso_name == "Association for Monitoring Equal Rights", 1, onlymostrecent$english_website)
onlymostrecent$english_website <- ifelse(onlymostrecent$cso_name == "Confederation of Public Employees Trade Unions (KESK)", 1, onlymostrecent$english_website)
onlymostrecent$english_website <- ifelse(onlymostrecent$cso_name == "Diyarbakir Institute for Political and Social Research", 1, onlymostrecent$english_website)
onlymostrecent$english_website <- ifelse(onlymostrecent$cso_name == "Turkish Medical association", 1, onlymostrecent$english_website)
onlymostrecent$english_website <- ifelse(onlymostrecent$cso_name == "Alliance for Finance Monitoring (ACFIM)", 1, onlymostrecent$english_website)
onlymostrecent$english_website <- ifelse(onlymostrecent$cso_name == "DefendDefenders (East and Horn of Africa Human Rights Defenders Project)", 1, onlymostrecent$english_website)

# Regression of time since first Communication on most recent rating
m1 <- lm(relation_gov_yearly_num ~ years_since, data = onlymostrecent)
summary(m1)
m2 <- lm(relation_gov_yearly_num ~ years_since + human_rights_cso + english_website, data = onlymostrecent)
summary(m2)
m3 <- lm(relation_gov_yearly_num ~ years_since + factor(cowc), data = onlymostrecent)
summary(m3)
m4 <- lm(relation_gov_yearly_num ~ years_since + human_rights_cso + english_website + factor(cowc), data = onlymostrecent)
summary(m4)

##### Table A.20 for Appendix D #####
star_full <- stargazer(m1, m2, m3, m4
                       , type="latex"
                       , title = "Reported CSO-government relations and UNSP communications"
                       , digits=3
                       , font.size="scriptsize"
                       , notes.align ="l"
                       , omit.stat = c("n", "ser", "rsq", "adj.rsq")
                       , add.lines = list(c("Number of observations"
                                            , length(m1$fitted.values)
                                            , length(m2$fitted.values)
                                            , length(m2$fitted.values)
                                            , length(m4$fitted.values)) ,
                                          c("Country fixed effects", "No", "No", "Yes", "Yes"),
                                          c("CSO-specific controls", "No", "Yes", "No", "Yes"))
                       , column.labels = c("(1)","(2)","(3)","(4)","(5)")
                       , dep.var.labels.include = F
                       , dep.var.caption = "DV: Self-reported relation with government (original survey)"
                       , column.sep.width = "1pt"
                       , align=F
                       , model.numbers=F, model.names=F
                       , label="tab:relation_gov_yearly_crosscsos"
                       , table.placement="H"
                       , star.char = custom_significance
                       , star.cutoffs = c(0.1, 0.05, 0.01, 0.001))


note.latex <- "\\multicolumn{6}{l} {\\parbox[t]{12cm}{ \\textit{Notes:} The models are estimated with ordinary least squares regression. Model 2 and 4 control for whether CSOs are primarily dedicated to human rights, and whether they have an English website. $^{\\dag}$p$<$0.1; $^{*}$p$<$0.05; $^{**}$p$<$0.01; $^{***}$p$<$0.001}.} \\\\"
star_full[grepl("Note",star_full)] <- note.latex
cat(star_full, sep = "\n")
#write(star_full, "./04_analysis_cross_csos/outputs/relation_gov_yearly_crosscso.tex")
write(star_full, "./10_replication_files/Outputs/Tables/TableA20_AppD_relation_gov_yearly_crosscso.tex")


# Regression of time since last Communication on most recent rating
m1 <- lm(relation_gov_yearly_num ~ since_last_mentioned, data = onlymostrecent)
summary(m1)
m2 <- lm(relation_gov_yearly_num ~ since_last_mentioned + human_rights_cso + english_website, data = onlymostrecent)
summary(m2)
m3 <- lm(relation_gov_yearly_num ~ since_last_mentioned + factor(cowc), data = onlymostrecent)
summary(m3)
m4 <- lm(relation_gov_yearly_num ~ since_last_mentioned + human_rights_cso + english_website + factor(cowc), data = onlymostrecent)
summary(m4)


star_full <- stargazer(m1, m2, m3, m4
                       , type="latex"
                       , title = "Reported CSO-government relations and UNSP communications"
                       , digits=3
                       , font.size="scriptsize"
                       , notes.align ="l"
                       , omit.stat = c("n", "ser", "rsq", "adj.rsq")
                       , add.lines = list(c("Number of observations"
                                            , length(m1$fitted.values)
                                            , length(m2$fitted.values)
                                            , length(m2$fitted.values)
                                            , length(m4$fitted.values)) ,
                                          c("Country fixed effects", "No", "No", "Yes", "Yes"),
                                          c("CSO-specific controls", "No", "Yes", "No", "Yes"))
                       , column.labels = c("(1)","(2)","(3)","(4)","(5)")
                       , dep.var.labels.include = F
                       , dep.var.caption = "DV: Self-reported relation with government (original survey)"
                       , column.sep.width = "1pt"
                       , align=F
                       , model.numbers=F, model.names=F
                       , label="tab:relation_gov_yearly_crosscsos"
                       , table.placement="H"
                       , star.char = custom_significance
                       , star.cutoffs = c(0.1, 0.05, 0.01, 0.001))


note.latex <- "\\multicolumn{6}{l} {\\parbox[t]{12cm}{ \\textit{Notes:} The models are estimated with ordinary least squares regression. Model 2 and 4 control for whether CSOs are primarily dedicated to human rights, and whether they have an English website. $^{\\dag}$p$<$0.1; $^{*}$p$<$0.05; $^{**}$p$<$0.01; $^{***}$p$<$0.001}.} \\\\"
star_full[grepl("Note",star_full)] <- note.latex
cat(star_full, sep = "\n")
#write(star_full, "./04_analysis_cross_csos/outputs/relation_gov_yearly_crosscso_lastmentioned.tex")


#### Subset analysis only with CSOs that remember whether a complaint was filed ####

# Create sub-dataset with only CSOs that remember whether a complaint was filed
ldata_remember <- ldata %>% filter(filed_known == "Yes")
controls <- paste0("+ v2peedueq_lag1 + conflict_lag + v2mecenefm_lag1 + v2juhcind_lag1 +  log_gdpcap_lag1")


# Binary variable for year of communication and all subsequent years
f0 <- as.formula( paste0("relation_gov_yearly_num  ~ treatment | 0 | 0 | cso_name + year"))
f1 <- as.formula( paste0("relation_gov_yearly_num ~  treatment | cso_name | 0 | cso_name + year"))
f2 <- as.formula( paste0("relation_gov_yearly_num ~  treatment ", controls,  " + years_since | cso_name | 0 | cso_name + year"))
f3 <- as.formula( paste0("relation_gov_yearly_num ~  treatment ", controls, " + years_since + lag(relation_gov_yearly_num,1) | lag(relation_gov_yearly_num,2:99)") )
f4 <- as.formula( paste0("relation_gov_yearly_num ~  " , controls, " + years_since + spvisits_lag1 | cso_name | (treatment ~  new_entries_sum + Rapporteurs_active_number) | cso_name + year"))

m0 <- felm(formula = f0, data = ldata_remember)
m1 <- felm(formula = f1, data = ldata_remember)
m2 <- felm(formula = f2, data = ldata_remember)
m3 <- pgmm(formula = f3, data = ldata_remember, index=c("cso_name", "year"), transformation = 'ld',  model = 'twostep', effect = 'individual')
m4 <- felm(formula = f4, data = ldata_remember)

# First stage
f4_fs <- as.formula(paste0("treatment ~ new_entries_sum + Rapporteurs_active_number", controls,  "+ spvisits_lag1 + years_since | cso_name | 0 | cso_name"))
m4_fs <- felm(formula = f4_fs, data = model.frame(m4))
summary(m4_fs)

# Summary stats
fstat <- round(summary(m4_fs)$fstat, 3)
fstat
sargan <- plm::sargan(m3)
sargan
autocorr1 <- mtest(m3,1)
autocorr1
autocorr2 <- mtest(m3,2)
autocorr2  
autocorr3 <- mtest(m3,3)
autocorr3  

##### Table A.19 for Appendix D #####
star_full <- stargazer(m0, m1, m2, m3, m4
                       , type="latex"
                       , title = "Reported CSO-government relations and UNSP communications"
                       , digits=3
                       , font.size="scriptsize"
                       , notes.align ="l"
                       , omit.stat = c("n", "ser", "rsq", "adj.rsq")
                       , add.lines = list(c("Number of observations"
                                            , length(m0$fitted.values)
                                            , length(m1$fitted.values)
                                            , length(m2$fitted.values)
                                            , length(m2$fitted.values)
                                            , length(m4$fitted.values)) ,
                                          c("Country fixed effects", "No", "Yes", "Yes", "Yes", "Yes"),
                                          c("Time-varying controls", "No", "No", "Yes", "Yes", "Yes"),
                                          c("Lagged depend. var.", "No", "No", "No", "Yes", "No"),
                                          c("First stage F-statistic", "NA", "NA", "NA", "NA",  fstat))
                       , covariate.labels = c("Post-communication years"
                                              , "Education (lag 1 yr)"
                                              , "Conflict (lag 1 yr)"
                                              , "Freedom of media (lag 1 yr)"
                                              , "Judicial independence (lag 1 yr)"
                                              , "GDP per capita (lag 1 yr)"
                                              , "Years since first communication"
                                              , "CSO-reported repression (lag 1 yr, instr.)"
                                              , "UNSP country visits"
                                              , "Post-communication (instr.)")
                       , column.labels = c("(1)","(2)","(3)","(4)","(5)")
                       , dep.var.labels.include = F
                       , dep.var.caption = "DV: Self-reported relation with government (original survey)"
                       , column.sep.width = "1pt"
                       , align=F
                       , model.numbers=F, model.names=F
                       , label="tab:relation_gov_yearly_full_subsetremember"
                       , table.placement="H"
                       , star.char = custom_significance
                       , star.cutoffs = c(0.1, 0.05, 0.01, 0.001))


note.latex <- "\\multicolumn{2}{l} {\\parbox[t]{14cm}{ \\textit{Notes:} The models in columns 1, 2, 3, and 5 are estimated with ordinary least squares regression, while the model in column 4 is estimated using generalized methods of moments. The lagged dependent variable in the model in column 4 is instrumented with deeper lags. Cluster-robust standard errors are in parentheses. $^{\\dag}$p$<$0.1; $^{*}$p$<$0.05; $^{**}$p$<$0.01; $^{***}$p$<$0.001}.} \\\\"
star_full[grepl("Note",star_full)] <- note.latex
cat(star_full, sep = "\n")
#write(star_full, "./04_analysis_cross_csos/outputs/relation_gov_yearly_subset_remember.tex")
write(star_full, "./10_replication_files/Outputs/Tables/TableA19_AppD_relation_gov_yearly_subset_remember.tex")


### exit #####

